#front matter
rm(list=ls())
library(foreign)
library(car)
library(lme4)
library(lattice)
library(arm)
library(texreg)

#set working directory
setwd('/Volumes/MONOGAN/immigration/demDef/')

#read data
data<-read.dta('healthImmigLONG.dta')
lp.data<-read.csv('LaxPhillipsAJPSmodified.csv')
lp.7<-subset(lp.data,subset=issue.full%in%c('bilingual','immlicense','tuition','verify','marijuanamedical','schipsixtyfull','suicide'))

###clean data###
#indicator for immigration issue
data$immig<-as.numeric(data$issue=="bilingual" | data$issue=="immlicense" | data$issue=="tuition" | data$issue=="verify")

#salience score
data$salience<-308+1372*data$immig

#create congruence measure
data$congruent<-as.numeric((data$liberalpolicy==1 & data$opinion>=50) | (data$liberalpolicy==0 & data$opinion<50))

#size of opinion majority
data$sizeMajority<-abs(data$opinion-50)

#conservative opinion majority
data$consMajority<-as.numeric(data$opinion<50)

#old variables to incorporate
data$old.scale.nyt<- -0.521374541
data$old.scale.nyt[data$issue%in%c('marijuanamedical','schip','suicide')]<- 0.03756527

#rescale continuous variables
data$Zsalience<-(data$salience-1092)/(2*679.9359)
data$Zpubideolunweighted<-rescale(data$pubideolunweighted)
data$Zgovt_ideol <-rescale(data$govt_ideol)
data$Zgov_avg<-rescale(data$gov_avg)
data$Zleg_avg <-rescale(data$leg_avg)
data$ZSquire2009 <-rescale(data$squire2009)
data$Zavg_turnout <-rescale(data$avg_turnout)
data$Zone_party_years <-rescale(data$one_party_years)
data$Zopinion<-rescale(data$opinion)
data$ZsizeMajority<-rescale(data$sizeMajority)
data$Zturnout12 <-rescale(data$turnout12)
data$Zone_party_years0914 <-rescale(data$one_party_years0914)
data$Zgovt_ideol0914 <-rescale(data$govt_ideol0914)
data$Zgov_avg0914<-rescale(data$gov_avg0914)
data$Zleg_avg0914 <-rescale(data$leg_avg0914)

#voter ideological opposition
data$voterIdeolOpp<-ifelse(data$opinion<50,data$Zpubideolunweighted,-data$Zpubideolunweighted) #If opinion is conservative, public liberalism. Else -1*liberalism.

#governmental ideological opposition
data$govIdeolOpp<-ifelse(data$opinion<50,data$Zgovt_ideol,-data$Zgovt_ideol) #if opinion is conservative, government liberalism. Else -1*liberalism.
data$govIdeolOpp0914<-ifelse(data$opinion<50,data$Zgovt_ideol0914,-data$Zgovt_ideol0914) #if opinion is conservative, government liberalism. Else -1*liberalism.

#legislative partisan opposition (L&P got this backwards)
data$legPartyOpp<-ifelse(data$opinion<50,data$Zleg_avg,-data$Zleg_avg) #if opinion is conservative, percent Democrat. Else -1*Democrat.
data$legPartyOpp0914<-ifelse(data$opinion<50,data$Zleg_avg0914,-data$Zleg_avg0914) #if opinion is conservative, percent Democrat. Else -1*Democrat.

#gubernatorial ideological opposition (L&P got this backwards)
data$govPartyOpp<-ifelse(data$opinion<50,data$Zgov_avg,-data$Zgov_avg) #if opinion is conservative, percent Democrat. Else -1*Democrat.
data$govPartyOpp0914<-ifelse(data$opinion<50,data$Zgov_avg0914,-data$Zgov_avg0914) #if opinion is conservative, percent Democrat. Else -1*Democrat.

#interest group opposition
data$igOp<-ifelse(data$opinion<50,data$powerfulig,-data$powerfulig) #if opinion is conservative, interest group liberalism. Else -1*liberalism.

###TABLE OF DEMOCRATIC DEFICIT BY POLICY###
#percent liberal policy
table(data$liberalpolicy[data$issue=="bilingual"])*2 #times 2 gives percentages for 50 states
table(data$liberalpolicy[data$issue=="immlicense"])*2
table(data$liberalpolicy[data$issue=="tuition"])*2
table(data$liberalpolicy[data$issue=="verify"])*2
table(data$liberalpolicy[data$issue=="marijuanamedical"])*2
table(data$liberalpolicy[data$issue=="schip"])*2
table(data$liberalpolicy[data$issue=="suicide"])*2

#average liberal opinion
round(mean(data$opinion[data$issue=="bilingual"]))
round(mean(data$opinion[data$issue=="immlicense"]))
round(mean(data$opinion[data$issue=="tuition"]))
round(mean(data$opinion[data$issue=="verify"]))
round(mean(data$opinion[data$issue=="marijuanamedical"]))
round(mean(data$opinion[data$issue=="schip"]))
round(mean(data$opinion[data$issue=="suicide"]))

#percent conservative opinion majorities (percent zeros=liberal, goes into Table)
table(data$consMajority[data$issue=="bilingual"])*2
table(data$consMajority[data$issue=="immlicense"])*2
table(data$consMajority[data$issue=="tuition"])*2
table(data$consMajority[data$issue=="verify"])*2
table(data$consMajority[data$issue=="marijuanamedical"])*2
table(data$consMajority[data$issue=="schip"])*2
table(data$consMajority[data$issue=="suicide"])*2

#percent congruence
table(data$congruent[data$issue=="bilingual"])*2
table(data$congruent[data$issue=="immlicense"])*2
table(data$congruent[data$issue=="tuition"])*2
table(data$congruent[data$issue=="verify"])*2
table(data$congruent[data$issue=="marijuanamedical"])*2
table(data$congruent[data$issue=="schip"])*2
table(data$congruent[data$issue=="suicide"])*2

#liberal incongruence bias, count of zeroes
table(data$congruent[data$liberalpolicy==1 &data$issue=="bilingual"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="immlicense"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="tuition"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="verify"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="marijuanamedical"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="schip"])
table(data$congruent[data$liberalpolicy==1 &data$issue=="suicide"])

#conservative incongruence bias, count of zeroes
table(data$congruent[data$liberalpolicy==0 &data$issue=="bilingual"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="immlicense"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="tuition"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="verify"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="marijuanamedical"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="schip"])
table(data$congruent[data$liberalpolicy==0 &data$issue=="suicide"])

#combined congruence rate
all.cong<-c(80,80,66,46,38,34,22);mean(all.cong) 

###GRAPHS OF PUBLIC SENTIMENT###
#setwd('/Volumes/MONOGAN/immigration/demDef/document/')

#sort data for graphs
data<-data[order(data$issue,data$state),]

#Drivers' Licenses
#pdf('oLicensesA.pdf',width=5,height=5,pointsize=12)
license<-cbind(data$opinion[data$issue=="immlicense"],recode(data$liberalpolicy[data$issue=="immlicense"],'0=1;1=17'))
rownames(license)<-data$abbr[1:50]
license<-license[order(license[,1]),]
dotplot(license[,1],pch=license[,2],col=as.factor(license[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Favoring Drivers' Licenses\n",scales=list(cex=.5))
#dev.off()

#E-Verify
#pdf('oEVerifyA.pdf',width=5,height=5,pointsize=12)
verify<-cbind(data$opinion[data$issue=="verify"],recode(data$liberalpolicy[data$issue=="verify"],'0=1;1=17'))
rownames(verify)<-data$abbr[1:50]
verify <-verify[order(verify[,1]),]
dotplot(verify[,1],pch= verify[,2],col=as.factor(verify[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Opposing E-Verify\n",scales=list(cex=.5))
#dev.off()

#In-State Tuition
#pdf('oTuitionA.pdf',width=5,height=5,pointsize=12)
tuition<-cbind(data$opinion[data$issue=="tuition"],recode(data$liberalpolicy[data$issue=="tuition"],'0=1;1=17'))
rownames(tuition)<-data$abbr[1:50]
tuition <-tuition[order(tuition[,1]),]
dotplot(tuition[,1],pch= tuition[,2],col=as.factor(tuition[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Favoring In-State Tuition\n",scales=list(cex=.5))
#dev.off()

#Bilingual Education
#pdf('oBilingualA.pdf',width=5,height=5,pointsize=12)
bilingual<-cbind(data$opinion[data$issue=="bilingual"],recode(data$liberalpolicy[data$issue=="bilingual"],'0=1;1=17'))
rownames(bilingual)<-data$abbr[1:50]
bilingual <-bilingual[order(bilingual[,1]),]
dotplot(bilingual[,1],pch= bilingual[,2],col=as.factor(bilingual[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Favoring Bilingual Education\n",scales=list(cex=.5))
#dev.off()

#Assisted Suicide
#pdf('oEuthanasiaA.pdf',width=5,height=5,pointsize=12)
suicide<-cbind(data$opinion[data$issue=="suicide"],recode(data$liberalpolicy[data$issue=="suicide"],'0=1;1=17'))
rownames(suicide)<-data$abbr[1:50]
suicide <-suicide[order(suicide[,1]),]
dotplot(suicide[,1],pch= suicide[,2],col=as.factor(suicide[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent to Allow Euthanasia\n",scales=list(cex=.5))
#dev.off()

#Medical Marijuana
#pdf('oMarijuanaA.pdf',width=5,height=5,pointsize=12)
marijuanamedical<-cbind(data$opinion[data$issue=="marijuanamedical"],recode(data$liberalpolicy[data$issue=="marijuanamedical"],'0=1;1=17'))
rownames(marijuanamedical)<-data$abbr[1:50]
marijuanamedical <-marijuanamedical[order(marijuanamedical[,1]),]
dotplot(marijuanamedical[,1],pch= marijuanamedical[,2],col=as.factor(marijuanamedical[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Favoring Medical Marijuana\n",scales=list(cex=.5))
#dev.off()

#SCHIP Coverage
#pdf('oSCHIPA.pdf',width=5,height=5,pointsize=12)
schip<-cbind(data$opinion[data$issue=="schip"],recode(data$liberalpolicy[data$issue=="schip"],'0=1;1=17'))
rownames(schip)<-data$abbr[1:50]
schip <-schip[order(schip[,1]),]
dotplot(schip[,1],pch= schip[,2],col=as.factor(schip[,2]),
key=list(text=list(c('Conservative Policy','Liberal Policy'),col=c(1,2)),points=list(pch=c(1,17),col=c(1,2)),space="bottom",columns=2),
xlab="Percent Favoring More SCHIP Coverage\n",scales=list(cex=.5))
#dev.off()

###ADD LAGS###
data$issue[data$issue=="schip"]<-"schipsixtyfull"
data<-data[order(data$issue,data$abbr),]
lp.7<-lp.7[order(lp.7$issue,lp.7$ostate.full),]
table(lp.7$issue==data$issue)
table(lp.7$ostate.full==data$abbr)
data$l.policy<-lp.7$have.liberal.policy
data$l.congruent<-lp.7$congruence

###DESCRIPTIVE STATISTICS###
summary(data)
summary(lp.data)

###PARTISANSHIP AND GOVERNMENT IDEOLOGY####
plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Bilingual Education")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="bilingual" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="bilingual" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="bilingual" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="bilingual" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Bilingual Education")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="bilingual" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="bilingual" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="bilingual" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="bilingual" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="Bilingual Education")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="bilingual" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="bilingual" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="bilingual" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="bilingual" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="Bilingual Education")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="bilingual" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="bilingual" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="bilingual" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="bilingual" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Immigrant Licenses")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="immlicense" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="immlicense" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="immlicense" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="immlicense" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Immigrant Licenses")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="immlicense" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="immlicense" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="immlicense" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="immlicense" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="Immigrant Licenses")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="immlicense" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="immlicense" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="immlicense" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="immlicense" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="Immigrant Licenses")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="immlicense" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="immlicense" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="immlicense" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="immlicense" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Medical Marijuana")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="marijuanamedical" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="marijuanamedical" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="marijuanamedical" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="marijuanamedical" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Medical Marijuana")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="marijuanamedical" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="marijuanamedical" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="marijuanamedical" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="marijuanamedical" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="Medical Marijuana")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="marijuanamedical" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="marijuanamedical" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="marijuanamedical" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="marijuanamedical" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="Medical Marijuana")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="marijuanamedical" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="marijuanamedical" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="marijuanamedical" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="marijuanamedical" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="SCHIP $60,0000")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="schipsixtyfull" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="schipsixtyfull" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="SCHIP $60,0000")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="schipsixtyfull" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="schipsixtyfull" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="SCHIP $60,0000")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="schipsixtyfull" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="schipsixtyfull" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="SCHIP $60,0000")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="schipsixtyfull" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="schipsixtyfull" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="schipsixtyfull" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Doctor-Assisted Suicide")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="suicide" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="suicide" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="suicide" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="suicide" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Doctor-Assisted Suicide")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="suicide" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="suicide" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="suicide" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="suicide" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="Doctor-Assisted Suicide")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="suicide" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="suicide" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="suicide" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="suicide" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="Doctor-Assisted Suicide")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="suicide" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="suicide" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="suicide" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="suicide" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Immigrant In-State Tuition")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="tuition" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="tuition" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="tuition" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="tuition" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="Immigrant In-State Tuition")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="tuition" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="tuition" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="tuition" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="tuition" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="Immigrant In-State Tuition")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="tuition" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="tuition" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="tuition" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="tuition" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="Immigrant In-State Tuition")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="tuition" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="tuition" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="tuition" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="tuition" & data$liberalpolicy==1],col='red',pch=3)


plot(y=data$Zgovt_ideol,x=data$Zgov_avg,type='n',xlab="Gubernatorial Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="E-Verify")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="verify" & data$liberalpolicy==0],x=data$Zgov_avg[data$issue=="verify" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="verify" & data$liberalpolicy==1],x=data$Zgov_avg[data$issue=="verify" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgovt_ideol,x=data$Zleg_avg,type='n',xlab="Legislative Control (Higher Is More Democratic)",ylab="Government Ideology (Higher is More Liberal)",main="E-Verify")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgovt_ideol[data$issue=="verify" & data$liberalpolicy==0],x=data$Zleg_avg[data$issue=="verify" & data$liberalpolicy==0])
points(y=data$Zgovt_ideol[data$issue=="verify" & data$liberalpolicy==1],x=data$Zleg_avg[data$issue=="verify" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zgov_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Gubernatorial Control (Higher Is More Democratic)",main="E-Verify")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zgov_avg[data$issue=="verify" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="verify" & data$liberalpolicy==0])
points(y=data$Zgov_avg[data$issue=="verify" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="verify" & data$liberalpolicy==1],col='red',pch=3)

plot(y=data$Zleg_avg,x=data$Zpubideolunweighted,type='n',xlab="Public Ideology (Higher is More Liberal)",ylab="Legislative Control (Higher Is More Democratic)",main="E-Verify")
abline(v=0,col='gray50');abline(h=0,col='gray50')
points(y=data$Zleg_avg[data$issue=="verify" & data$liberalpolicy==0],x=data$Zpubideolunweighted[data$issue=="verify" & data$liberalpolicy==0])
points(y=data$Zleg_avg[data$issue=="verify" & data$liberalpolicy==1],x=data$Zpubideolunweighted[data$issue=="verify" & data$liberalpolicy==1],col='red',pch=3)




###CORRELATION BETWEEN GOVERNMENT IDEOLOGY AND ELITE PARTISANSHIP###
cor(lp.7$rescale.Berry.full.,lp.7$rescale.legislative.share.dem.mean.full.,use="pairwise")
cor(lp.7$rescale.Berry.full.,lp.7$rescale.gov.control.mean.dem.full.,use="pairwise")
cor(data$Zgovt_ideol,data$Zleg_avg,use="pairwise")
cor(data$Zgovt_ideol,data$Zgov_avg,use="pairwise")



###RESPONSIVENESS MODELS###
#Lax & Phillips 2012, Table 3, Model 6: pure replication
lp.full.table3.mod6.pure<-glmer(have.liberal.policy~rescale.opinion.*rescale.nyt.bypolicygroup.full.+rescale.EWM.full.+rescale.Berry.full.+rescale.legislative.share.dem.mean.full.+rescale.gov.control.mean.dem.full.+rescale.prof.full.*rescale.opinion.+term.limits.full*rescale.opinion.+dirdem.full*rescale.opinion.+courtelected.full*rescale.opinion.+Powerful.full.strong.balance+rescale.turnout.full.*rescale.opinion.+rescale.onepartyness.full.*rescale.opinion.+(rescale.opinion.|ostate.full)+(rescale.opinion.|issue.full),family=binomial(link='logit'),data=lp.data);summary(lp.full.table3.mod6.pure)
p.1<-as.numeric(predict(lp.full.table3.mod6.pure,type="response")>0.5)
mean(p.1==lp.data$have.liberal.policy[lp.data$ostate.full!="NE"])
table(p.1,lp.data$have.liberal.policy[lp.data$ostate.full!="NE"])

#Same model with Lax & Phillips data subsetted to seven issues under study
lp.7issue.table3.mod6<-glmer(have.liberal.policy~rescale.opinion.*rescale.nyt.bypolicygroup.full.+rescale.EWM.full.+rescale.Berry.full.+rescale.legislative.share.dem.mean.full.+rescale.gov.control.mean.dem.full.+rescale.prof.full.*rescale.opinion.+term.limits.full*rescale.opinion.+dirdem.full*rescale.opinion.+courtelected.full*rescale.opinion.+Powerful.full.strong.balance+rescale.turnout.full.*rescale.opinion.+rescale.onepartyness.full.*rescale.opinion.+(rescale.opinion.|ostate.full)+(rescale.opinion.|issue.full),family=binomial(link='logit'),data=lp.7);summary(lp.7issue.table3.mod6)
p.2<-as.numeric(predict(lp.7issue.table3.mod6,type="response")>0.5)
mean(p.2==lp.7$have.liberal.policy[lp.7$ostate.full!="NE"])
table(p.2,lp.7$have.liberal.policy[lp.7$ostate.full!="NE"])

#Model with updated data
new.table3.mod6<-glmer(liberalpolicy~ Zopinion*Zsalience+Zpubideolunweighted+Zgovt_ideol+Zleg_avg+Zgov_avg+ZSquire2009*Zopinion+termlimit*Zopinion+initiative*Zopinion+electedcourt*Zopinion+powerfulig+Zavg_turnout*Zopinion+Zone_party_years*Zopinion+(Zopinion|abbr)+(Zopinion|issue),family=binomial(link='logit'),data=data,subset=nebraska==0);summary(new.table3.mod6)
p.3<-as.numeric(predict(new.table3.mod6,type="response")>0.5)
mean(p.3==data$liberalpolicy[data$nebraska==0])
table(p.3,data$liberalpolicy[data$nebraska==0])

#Model with updated data and lag
new.table3.mod6.lag<-glmer(liberalpolicy~ l.policy+Zopinion*Zsalience+Zpubideolunweighted+Zgovt_ideol+Zleg_avg+Zgov_avg+ZSquire2009*Zopinion+termlimit*Zopinion+initiative*Zopinion+electedcourt*Zopinion+powerfulig+Zavg_turnout*Zopinion+Zone_party_years*Zopinion+(Zopinion|abbr)+(Zopinion|issue),family=binomial(link='logit'),data=data,subset=nebraska==0);summary(new.table3.mod6.lag)
p.4<-as.numeric(predict(new.table3.mod6.lag,type="response")>0.5)
mean(p.4==data$liberalpolicy[data$nebraska==0])
table(p.4,data$liberalpolicy[data$nebraska==0])


##output to LaTeX##
texreg(list(lp.7issue.table3.mod6,new.table3.mod6,new.table3.mod6.lag),custom.model.names=c('2008','2014','2014 with Lag'),custom.coef.names=c('Intercept','Opinion','Salience','Voter Liberalism','Govt. Liberalism','Democratic Legislature Perc.','Democratic Governor','Professionalization','Term Limits','Citizen Init.','Elected Court','Powerful Int. Group Balance','Turnout','One-Party Dominance','Salience$\\times$Op.','Professionalization$\\times$Op.','Term Limits$\\times$Op.','Citizen Init.$\\times$Op.','Elected Court$\\times$Op.','Turnout$\\times$Op.','One-Party Dominance$\\times$Op.',  'Opinion', 'Salience', 'Voter Liberalism','Govt. Liberalism','Democratic Legislature Perc.','Democratic Governor','Professionalization','Term Limits','Citizen Init.','Elected Court','Powerful Int. Group Balance','Turnout','One-Party Dominance','Salience$\\times$Op.','Professionalization$\\times$Op.','Term Limits$\\times$Op.','Citizen Init.$\\times$Op.','Elected Court$\\times$Op.','Turnout$\\times$Op.','One-Party Dominance$\\times$Op.','Lagged policy'),stars=c(.1),dcolumn=T,label="responsiveness",caption="\\textbf{Policy Responsiveness Models}",caption.above=T,custom.gof.names=c('AIC','BIC','Log Likelihood','Num.~obs.','Num.~states','Num.~issues','Variance: state intercepts','Variance: state opinion slopes','Covariance in state rand.~eff.','Variance: issue intercepts','Variance: issue opinion slopes','Covariance in issue rand.~eff.','Num.~states','Num.~issues','Variance: state intercepts','Variance: state opinion slopes','Covariance in state rand.~eff.','Variance: issue intercepts','Variance: issue opinion slopes','Covariance in issue rand.~eff.'))

texreg(list(lp.full.table3.mod6.pure),custom.model.names=c('39 Issues 2008'),custom.coef.names=c('Intercept','Opinion','Salience','Voter Liberalism','Govt. Liberalism','Democratic Legislature Perc.','Democratic Governor','Professionalization','Term Limits','Citizen Init.','Elected Court','Powerful Int. Group Balance','Turnout','One-Party Dominance','Salience$\\times$Op.','Professionalization$\\times$Op.','Term Limits$\\times$Op.','Citizen Init.$\\times$Op.','Elected Court$\\times$Op.','Turnout$\\times$Op.','One-Party Dominance$\\times$Op.'),stars=c(.1),dcolumn=T,label="responsivenessLP",caption="\\textbf{Policy Responsiveness Models}",caption.above=T)

# Run interactions one-at-a-time
source("interactions.oneatatime.R",echo=T)

###CONGRUENCE MODELS###
#Lax & Phillips 2012, Table 4, Model 6: pure replication
lp.full.table4.mod6.pure<-glmer(congruence~rescale.size.of.majority.+rescale.nyt.bypolicygroup.full.+conservative.opinion.majority+flip.EWM.against+flip.Berry.against+flip.legislative.share.dem.mean.full+flip.gov.control.mean.dem.full+rescale.prof.full.+term.limits.full+dirdem.full+courtelected.full+flip.powerful.balance.against+rescale.turnout.full.+rescale.onepartyness.full.+(rescale.size.of.majority.|ostate.full)+(rescale.size.of.majority.|issue.full),family=binomial(link='logit'),data=lp.data);summary(lp.full.table4.mod6.pure)
q.1<-as.numeric(predict(lp.full.table4.mod6.pure,type="response")>0.5)
mean(q.1==lp.data$congruence[lp.data$ostate.full!="NE"])
table(q.1,lp.data$congruence[lp.data$ostate.full!="NE"])

#Same model with Lax & Phillips data subsetted to seven issues under study
lp.7issue.table4.mod6<-glmer(congruence~rescale.size.of.majority.+rescale.nyt.bypolicygroup.full.+conservative.opinion.majority+flip.EWM.against+flip.Berry.against+flip.legislative.share.dem.mean.full+flip.gov.control.mean.dem.full+rescale.prof.full.+term.limits.full+dirdem.full+courtelected.full+flip.powerful.balance.against+rescale.turnout.full.+rescale.onepartyness.full.+(rescale.size.of.majority.|ostate.full)+(rescale.size.of.majority.|issue.full),family=binomial(link='logit'),data=lp.7);summary(lp.7issue.table4.mod6)
q.2<-as.numeric(predict(lp.7issue.table4.mod6,type="response")>0.5)
mean(q.2==lp.7$congruence[lp.7$ostate.full!="NE"])
table(q.2,lp.7$congruence[lp.7$ostate.full!="NE"])

#Model with updated data
new.table4.mod6<-glmer(congruent~ZsizeMajority+Zsalience+consMajority+voterIdeolOpp+govIdeolOpp+legPartyOpp+govPartyOpp+ZSquire2009+termlimit+initiative+electedcourt+igOp+Zavg_turnout+Zone_party_years+(ZsizeMajority|abbr)+(ZsizeMajority|issue),family=binomial(link='logit'),data=data,subset=nebraska==0);summary(new.table4.mod6)
q.3<-as.numeric(predict(new.table4.mod6,type="response")>0.5)
mean(q.3==data$congruent[data$nebraska==0])
table(q.3,data$congruent[data$nebraska==0])

#Model with updated data and lag
new.table4.mod6.lag<-glmer(congruent~l.congruent+ZsizeMajority+Zsalience+consMajority+voterIdeolOpp+govIdeolOpp+legPartyOpp+govPartyOpp+ZSquire2009+termlimit+initiative+electedcourt+igOp+Zavg_turnout+Zone_party_years+(ZsizeMajority|abbr)+(ZsizeMajority|issue),family=binomial(link='logit'),data=data,subset=nebraska==0);summary(new.table4.mod6.lag)
q.4<-as.numeric(predict(new.table4.mod6.lag,type="response")>0.5)
mean(q.4==data$congruent[data$nebraska==0])
table(q.4,data$congruent[data$nebraska==0])

#BLUP for two static models#
re.08<-ranef(lp.7issue.table4.mod6)
re.14<-ranef(new.table4.mod6)

d.08<-data.frame(y08=re.08$ostate.full$"rescale.size.of.majority",abbr=rownames(re.08$ostate.full))
d.14<-data.frame(y14=re.14$abbr$"ZsizeMajority",abbr=rownames(re.14$abbr))

re.data<-merge(x=d.08,y=d.14,by="abbr")
re.data<-re.data[order(re.data$y14),]

pdf('re2014.pdf',width=5,height=5,pointsize=12)
plot(x=re.data$y14,y=1:49,type='n',xlab="State-Specific Random Slope Term 2014",ylab="",axes=F,xlim=c(-.58,.66))
axis(1)
axis(2,at=1:49,labels=re.data$abbr,las=2,cex.axis=.5)
abline(h=1:49,col='gray90')
box()
points(x=re.data$y14,y=1:49,pch=3,col='red')
dev.off()

pdf('re2008.pdf',width=5,height=5,pointsize=12)
options(scipen=12)
plot(x=re.data$y08,y=1:49,type='n',xlab="State-Specific Random Slope Term 2008",ylab="",axes=F,xlim=c(-.58,.66))
axis(1)
axis(2,at=1:49,labels=re.data$abbr,las=2,cex.axis=.5)
abline(h=1:49,col='gray90')
box()
points(x=re.data$y08,y=1:49,pch=20)
dev.off()

##output to LaTeX##
texreg(list(lp.7issue.table4.mod6,new.table4.mod6,new.table4.mod6.lag),custom.model.names=c('2008','2014','2014 with Lag'),custom.coef.names=c('Intercept','Size of Opinion Majority','Salience','Conservative Opinion Majority','Voter Ideological Opposition','Govt. Ideological Opposition','Legislative Partisan Opposition','Governor Partisan Opposition','Professionalization','Term Limits','Citizen Init.','Elected Court','Interest Group Opposition','Turnout','One-Party Dominance','Size of Opinion Majority','Salience','Conservative Opinion Majority','Voter Ideological Opposition','Govt. Ideological Opposition','Legislative Partisan Opposition','Governor Partisan Opposition','Professionalization','Term Limits','Citizen Init.','Elected Court','Interest Group Opposition','Turnout','One-Party Dominance','Lagged Congruence'),stars=c(.1),dcolumn=T,label="congruence",caption="\\textbf{Policy Congruence Models}",caption.above=T,custom.gof.names=c('AIC','BIC','Log Likelihood','Num.~obs.','Num.~states','Num.~issues','Variance: state intercepts','Variance: state opinion slopes','Covariance in state rand.~eff.','Variance: issue intercepts','Variance: issue opinion slopes','Covariance in issue rand.~eff.','Num.~states','Num.~issues','Variance: state intercepts','Variance: state opinion slopes','Covariance in state rand.~eff.','Variance: issue intercepts','Variance: issue opinion slopes','Covariance in issue rand.~eff.'))


texreg(list(lp.full.table4.mod6.pure),custom.model.names=c('39 Issues 2008'),custom.coef.names=c('Intercept','Size of Opinion Majority','Salience','Conservative Opinion Majority','Voter Ideological Opposition','Govt. Ideological Opposition','Legislative Partisan Opposition','Governor Partisan Opposition','Professionalization','Term Limits','Citizen Init.','Elected Court','Interest Group Opposition','Turnout','One-Party Dominance'),stars=c(.1),dcolumn=T,label="congruenceLP",caption="\\textbf{Policy Congruence Models}",caption.above=T)

###difference of proportion test###
#p = (p1 * n1 + p2 * n2) / (n1 + n2)
#SE = sqrt{ p * ( 1 - p ) * [ (1/n1) + (1/n2) ] }
#z = (p1 - p2) / SE
p<-(.52*350+.42*350)/(700);p
SE<-sqrt(.47*.53*((1/350)+(1/350)));SE
.1/SE
1-pnorm(.1/SE)


###ALTERNATIVE VERSIONS OF MODELS THAT USE ONLY POST-2008 INFORMATION###
#responsiveness
responsive.recent<-glmer(liberalpolicy~ l.policy+Zopinion*Zsalience+Zpubideolunweighted+Zgovt_ideol0914+Zleg_avg0914+Zgov_avg0914+ZSquire2009*Zopinion+termlimit*Zopinion+initiative*Zopinion+electedcourt*Zopinion+powerfulig+Zturnout12*Zopinion+Zone_party_years0914*Zopinion+(Zopinion|abbr)+(Zopinion|issue),family=binomial(link='logit'),data=data,subset=nebraska==0);summary(responsive.recent)
p.5<-as.numeric(predict(responsive.recent,type="response")>0.5)
mean(p.5==data$liberalpolicy[data$nebraska==0])
table(p.5,data$liberalpolicy[data$nebraska==0])

texreg(list(responsive.recent),custom.model.names=c('2014 Responsiveness'),custom.coef.names=c('Intercept','Lagged policy','Opinion','Salience','Voter Liberalism','Govt. Liberalism','Democratic Legislature Perc.','Democratic Governor','Professionalization','Term Limits','Citizen Init.','Elected Court','Powerful Int. Group Balance','Turnout','One-Party Dominance','Salience$\\times$Op.','Professionalization$\\times$Op.','Term Limits$\\times$Op.','Citizen Init.$\\times$Op.','Elected Court$\\times$Op.','Turnout$\\times$Op.','One-Party Dominance$\\times$Op.'),stars=c(.1),dcolumn=T,label="responsivenessRecent",caption="\\textbf{Policy Responsiveness Model for 2014 with Measures Computed 2009-2014}",caption.above=T)

#congruence
congruent.recent<-glmer(congruent~l.congruent+ZsizeMajority+Zsalience+consMajority+voterIdeolOpp+govIdeolOpp0914+legPartyOpp0914+govPartyOpp0914+ZSquire2009+termlimit+initiative+electedcourt+igOp+Zturnout12+Zone_party_years0914+(ZsizeMajority|abbr)+(ZsizeMajority|issue),family=binomial(link='logit'),data=data,subset=nebraska==0,nAGQ=0);summary(congruent.recent)#Fixed by: abbr slope, issue slope, dropping govt ideol, leg party, gov party, one party years, turnout, OR nAGQ=0 which is the option used. Still busted: lag
q.5<-as.numeric(predict(congruent.recent,type="response")>0.5)
mean(q.5==data$congruent[data$nebraska==0])
table(q.5,data$congruent[data$nebraska==0])

texreg(list(congruent.recent),custom.model.names=c('2014 Congruence'),custom.coef.names=c('Intercept','Lagged Congruence','Size of Opinion Majority','Salience','Conservative Opinion Majority','Voter Ideological Opposition','Govt. Ideological Opposition','Legislative Partisan Opposition','Governor Partisan Opposition','Professionalization','Term Limits','Citizen Init.','Elected Court','Interest Group Opposition','Turnout','One-Party Dominance'),stars=c(.1),dcolumn=T,label="congruenceRecent",caption="\\textbf{Policy Congruence Model for 2014 with Measures Computed 2009-2014}",caption.above=T)

